library(tidyverse)
library(readr)
library(readxl)
library(janitor)
library(stringr)
library(scales)
library(lmtest)
library(sandwich)
library(broom)
library(tibble)
library(jsonlite)
library(forcats)
library(ggplot2)
library(ggrepel)
library(gganimate)
library(EEAaq)
library(gifski)

1 Introducción

La fertilidad humana es sensible a condiciones ambientales que afectan tanto a la salud reproductiva como a las decisiones familiares.

En particular, la calidad del aire se ha relacionado con alteraciones hormonales, inflamación sistémica y estrés oxidativo que pueden reducir la fecundidad, aumentar el riesgo de pérdida gestacional y acortar la ventana fértil. Más allá de los mecanismos biológicos, el entorno en que las familias planifican también responde a expectativas de salud pública, seguridad ambiental y confianza institucional.

En este contexto, la vigilancia de la calidad del aire (monitorización sistemática, cobertura de estaciones, reporte y disponibilidad de datos) puede operar como un bien público que mejora la salud y reduce la incertidumbre sobre riesgos ambientales. Una red de vigilancia más densa y activa suele acompañar estrategias de control, cumplimiento normativo y comunicación del riesgo, y puede facilitar políticas de reducción de emisiones.

El resultado esperado (si estos canales funcionan) es un entorno más saludable y predecible para la crianza, lo que, en teoría, podría favorecer niveles más altos de fertilidad.

Sin embargo, a nivel agregado (país), se conoce menos acerca de:

Hasta qué punto la intensidad de la vigilancia del aire se asocia con niveles de fertilidad.
Si la contaminación atmosférica (por ejemplo, con NO₂ o CO₂) se relaciona con la fertilidad media de los países europeos.

El presente seminario aborda estas cuestiones desde una perspectiva descriptiva, utilizando datos agregados por país para Europa.

2 Objetivo general

Analizar la vigilancia y medición de la calidad del aire en relación con la fertilidad media (TFR) en los países europeos.

3 Objetivos específicos

Vigilancia vs fertilidad
¿Existe asociación entre la intensidad de la vigilancia de la calidad del aire y la fertilidad media (TFR) de los países europeos?

NO₂ vs fertilidad
¿Están los niveles medios de NO₂ en los países europeos relacionados con sus niveles medios de fertilidad (TFR)? Si es así, ¿de dónde proviene dicho NO₂?

CO₂ vs fertilidad
¿Existe una relación directamente proporcional entre los niveles de dióxido de carbono (CO₂) y la fertilidad?

PM₂.₅ en la vigilancia vs fertilidad
¿Se asocia la proporción de observaciones dedicadas a PM₂.₅ en los sistemas de vigilancia de la calidad del aire de los países europeos con su fertilidad media (TFR)?

4 Metodología y resultados:

4.1 Describir la vigilancia de la calidad del aire en Europa y ver su asociación con la fertilidad

Pregunta: ¿Existe una asociación entre la intensidad de vigilancia de la calidad del aire y la fertilidad en Europa?

Hipótesis: Se espera que una mayor vigilancia del aire implique una mayor fertilidad. Por ello, esperamos que los países con más observaciones de calidad del aire tengan una TFR ligeramente superior, aunque esta relación no implica causalidad directa; es un test descriptivo de la asociación.

Si la hipótesis es consistente con los datos, esperamos observar una relación positiva entre el nivel de vigilancia (más observaciones) y la TFR.

Nota: La vigilancia se mide como el número de observaciones registradas por país y año en el sistema de AQ e-Reporting (EEA). La fertilidad se mide mediante la TFR media anual por país (Eurostat).

4.1.1 Lectura y limpieza de Datos Base

Trabajamos con dos fuentes principales de datos:

Calidad del aire: obtenida mediante EEAaq::aq_annual_data() o CSV original (DataExtract.csv)
Fertilidad: obtenida de Eurostat (Fertilidad.csv)

Se leen los archivos, se realiza una limpieza inicial de columnas y se transforma el formato de los datos de Eurostat (códigos geo) al formato de nombres de país homogéneos (pais) para permitir la unión con los datos de vigilancia.

# 1. Lectura de datos y limpieza básica

aire <- read_delim(
  "INPUT/DATA/DataExtract.csv",
  delim = ",", escape_double = FALSE, trim_ws = TRUE
)

fertilidad <- read_delim(
  "INPUT/DATA/Fertilidad.csv",
  delim = ",", escape_double = FALSE, trim_ws = TRUE
)

aire       <- aire %>% clean_names()
fertilidad <- fertilidad %>% clean_names()

# 2. Tabla de mapeo de códigos de país (Eurostat -> nombre)

mapa_paises_df <- tibble::tribble(
  ~geo,        ~pais,
  "AD", "Andorra",
  "AL", "Albania",
  "AM", "Armenia",
  "AT", "Austria",
  "AZ", "Azerbaijan",
  "BA", "Bosnia and Herzegovina",
  "BE", "Belgium",
  "BG", "Bulgaria",
  "BY", "Belarus",
  "CH", "Switzerland",
  "CY", "Cyprus",
  "CZ", "Czechia",
  "DE", "Germany",
  "DK", "Denmark",
  "EE", "Estonia",
  "EL", "Greece",
  "ES", "Spain",
  "FI", "Finland",
  "FR", "France",
  "GE", "Georgia",
  "GI", "Gibraltar",
  "HR", "Croatia",
  "HU", "Hungary",
  "IE", "Ireland",
  "IS", "Iceland",
  "IT", "Italy",
  "LI", "Liechtenstein",
  "LT", "Lithuania",
  "LU", "Luxembourg",
  "LV", "Latvia",
  "MC", "Monaco",
  "MD", "Moldova",
  "ME", "Montenegro",
  "MK", "North Macedonia",
  "MT", "Malta",
  "NL", "Netherlands",
  "NO", "Norway",
  "PL", "Poland",
  "PT", "Portugal",
  "RO", "Romania",
  "RS", "Serbia",
  "RU", "Russia",
  "SE", "Sweden",
  "SI", "Slovenia",
  "SK", "Slovakia",
  "SM", "San Marino",
  "TR", "Türkiye",
  "UA", "Ukraine",
  "UK", "United Kingdom",
  "VA", "Vatican City",
  "XK", "Kosovo"
  # Agregados tipo "EU27_2020", etc. no se incluyen
)

4.1.2 Homogeneización de países y resumen de datos de fertilidad y aire

Se preparan los datos para los análisis posteriores unificando los nombres de los países entre distintas fuentes, de manera que los datasets de fertilidad y calidad del aire sean comparables. Se realizan varias tareas:

Homogeneización de nombres de países: Se recodifican alias o variantes de nombres para asegurar consistencia entre datasets.
Resumen de fertilidad: Se calcula la TFR media y el número de años con datos disponibles para cada país.
Conteo de observaciones de aire: Se contabilizan las observaciones de PM2.5 por país y año, aplicando la homogeneización de nombres.
Unión de datos: Se combinan los resúmenes de fertilidad y de observaciones de aire por país, obteniendo un dataset base que servirá para posteriores normalizaciones y análisis.

Este paso es clave para que los análisis estadísticos posteriores puedan comparar correctamente la vigilancia de la calidad del aire con la fertilidad a nivel país, evitando inconsistencias debidas a diferencias de nombres o códigos de país.

# Homogeneizar alias de países
recode_pais <- function(x) {
  x <- str_trim(x)
  dplyr::recode(
    x,
    "Turkey"              = "Türkiye",
    "Czech Republic"      = "Czechia",
    "Russian Federation"  = "Russia",
    "Republic of Moldova" = "Moldova",
    "UK"                  = "United Kingdom",
    "Great Britain"       = "United Kingdom",
    .default = x
  )
}

# Calcular la TFR media y número de años observados por país
resumir_fertilidad <- function(df) {
  df %>%
    group_by(pais) %>%
    summarise(
      n_anios    = n(),
      tfr_media  = mean(tfr, na.rm = TRUE),
      .groups    = "drop"
    ) %>%
    arrange(pais)
}

# Conteo de observaciones de aire por país-año y homogeneización de nombres
limpiar_y_contar_aire <- function(aire_raw) {
  aire_raw %>%
    mutate(
      pais = str_trim(country),
      anio = as.integer(year),
      # Aplicar homogeneización
      pais = recode_pais(country)) %>% 
    count(pais, anio, name = "n_obs_aire") %>%
    arrange(pais, anio)
}

# Filtrar y transformar los datos de Eurostat de código a nombre de país
fert_codigos <- fertilidad %>%
  filter(freq == "A") %>%
  transmute(geo = geo, anio = as.integer(time_period), tfr = as.numeric(obs_value))

fert_paises <- fert_codigos %>%
  left_join(mapa_paises_df, by = "geo") %>%
  filter(!is.na(pais)) %>%
  select(pais, anio, tfr)

# Resumen: TFR Media por país
fert_resumen_nombres <- resumir_fertilidad(fert_paises)

# Conteo de observaciones por país-año (Frecuencia de Vigilancia)
aire_conteo <- limpiar_y_contar_aire(aire)

# Resumen: Total de Observaciones por país (para normalización)
aire_resumen_pais <- aire_conteo %>%
  group_by(pais) %>%
  summarise(
    total_obs_aire      = sum(n_obs_aire, na.rm = TRUE),
    .groups             = "drop"
  )

# Resumen base combinado (para posterior normalización)
resumen_aire_fert_pais <- aire_resumen_pais %>%
  full_join(fert_resumen_nombres, by = "pais") %>%
  arrange(pais)

4.1.3 Normalización y Construcción de la Variable de Vigilancia

Para obtener una medida de vigilancia comparable entre países, es crucial normalizar por el tamaño poblacional. El proceso se realiza en dos fases:

Integración de la Infraestructura de Vigilancia: Se accede a los metadatos de las estaciones de la Agencia Europea de Medio Ambiente (EEA) a través del paquete EEAaq, obteniendo el número de estaciones activas por país (n_estaciones_activas), que es una medida de la infraestructura de monitorización.
Cálculo de Métricas Normalizadas: Se unen los datos de población (Población_Eurostat_GIND) a las métricas de vigilancia. Finalmente, se crean las dos variables principales:

Estaciones por millón de habitantes (estaciones_por_millon).
Observaciones por millón de habitantes (Log-transformadas) (obs_por_millon_log), siendo esta última la variable predictora clave (X) en el modelo de regresión.

El dataframe resultante, df_vigilancia_tfr, contendrá las métricas robustas de TFR y Vigilancia, listas para el análisis estadístico y la visualización.

# Obtener metadatos de estaciones de la EEA (dataframe interno de EEAaq)
vigilancia_infraestructura <- EEAaq_get_dataframe("stations") %>%
  clean_names() %>%
  mutate(
    pais = recode_pais(country)
  ) %>%
  # Contar estaciones únicas (infraestructura) por país
  group_by(pais) %>%
  summarise(
    n_estaciones_activas = n_distinct(air_quality_station_eo_i_code),
    .groups = "drop"
  )

# Carga de Datos de Población Reales (Usando el CSV de Eurostat DEMO_GIND) 

poblacion_raw <- read_delim(
  "INPUT/DATA/Poblacion_Eurostat_GIND.csv",
  delim = ",", escape_double = FALSE, trim_ws = TRUE
) %>%
  clean_names()

# Transformación de los datos de población:
poblacion_data <- poblacion_raw %>%
  
  # Filtramos por el año central de nuestro análisis
  filter(time_period == 2020) %>% 
  
  # Seleccionamos y transformamos las columnas clave
  transmute(
    country_name_raw = geo,
    poblacion_millones = as.numeric(obs_value) / 1000000 
  ) %>%
  
  # Homogeneizamos el nombre del país para el join
  mutate(pais = recode_pais(country_name_raw)) %>% 
  
  select(pais, poblacion_millones) %>%
  filter(!is.na(poblacion_millones))

# Unión y Cálculo de Métricas Normalizadas 

df_vigilancia_tfr <- resumen_aire_fert_pais %>%
  left_join(vigilancia_infraestructura, by = "pais") %>%
  left_join(poblacion_data, by = "pais") %>%
  
  # Eliminar filas donde el denominador o numerador sea NA/cero
  filter(
      poblacion_millones > 0,
      !is.na(total_obs_aire)
  ) %>%
  
  mutate(
    estaciones_por_millon = n_estaciones_activas / poblacion_millones,
    
    # Cálculo de la variable clave: la función log1p maneja el log(1+0) = 0 sin error.
    obs_por_millon_log = log1p(total_obs_aire / poblacion_millones)
  ) %>%
  
  # Para solo conservar filas con valores finitos
  filter(is.finite(tfr_media), is.finite(obs_por_millon_log))

4.1.4 Análisis Descriptivo I: Gráfico de barras

Se genera un gráfico de barras que visualiza la métrica de Infraestructura de Vigilancia normalizada (estaciones_por_millon). Para facilitar la visualización de la hipótesis, los países se ordenan por su TFR media, y el color de las barras se mapea a la TFR. Esto permite evaluar si una mayor TFR se asocia con una mayor altura de barra.

df_ordenado <- df_vigilancia_tfr %>%
  # Se crea un factor para ordenar el eje X por TFR_media ascendente
  mutate(pais_ordenado = forcats::fct_reorder(pais, tfr_media))
  
ggplot(df_ordenado, 
       aes(x = pais_ordenado, y = estaciones_por_millon)) +
  
  # Barras (geom_col), coloreadas por la TFR media
  geom_col(aes(fill = tfr_media), stat = "identity") +
  
  # Escala de color para representar la fertilidad
  scale_fill_gradient(
    low = "peachpuff", 
    high = "red", 
    space = "Lab", 
    name = "TFR media"
  ) +
  
  labs(
    title = "Asociación entre Infraestructura de Vigilancia (EEAaq) y Fertilidad",
    subtitle = "Países ordenados por TFR media. Altura = Estaciones de monitorización por millón de habitantes.",
    x = NULL, 
    y = "Nº de Estaciones Activas (por millón de hab.)"
  ) +
  
  scale_y_continuous(expand = expansion(mult = c(0, .1)), 
                     labels = scales::comma) +
  theme_classic() +
  theme(
    axis.text.x = element_text(angle = 60, vjust = 1, hjust = 1, size = 8),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

Resultados del gráfico de barras

El gráfico está ordenado de izquierda a derecha por la fertilidad (TFR) de los países. Los países con la TFR más baja están a la izquierda y los que tienen la TFR más alta están a la derecha.

Si existiera una asociación positiva (mayor vigilancia = mayor fertilidad), esperaríamos ver un patrón en el que las barras (altura) se hicieran progresivamente más altas a medida que se avanza de izquierda a derecha.

En la gráfica observamos que no hay un patrón claro en la altura de las barras a lo largo del eje x, lo que significa que no existe una relación entre la fertilidad con la mayor vigilancia.

Ejemplo: Países como Lituania y Serbia muestran barras extremadamente altas (muy alta vigilancia por habitante), pero su TFR (color) no es la más alta de Europa.

Refutación de la Hipótesis Inicial: El gráfico sugiere que la hipótesis de que “una mayor inversión en vigilancia se asocia con niveles de fertilidad ligeramente superiores” es probablemente falsa.

4.1.5 Análisis Descriptivo II: Gráfico de Frecuencia (Dispersión)

El diagrama de dispersión analiza la relación entre la Frecuencia de Reporte normalizada (obs_por_millon_log) y la TFR. Se incluye una recta de regresión lineal (OLS) para evaluar visualmente la tendencia y bandas de confianza para la incertidumbre.

df_etiquetas <- df_vigilancia_tfr %>%
  # Etiquetar países con valores extremos para control visual de outliers
  mutate(
    Etiqueta_pais = ifelse(
      tfr_media < 1.35 | tfr_media > 1.8 | obs_por_millon_log > 10,
      pais,
      ""
    )
  )

ggplot(df_etiquetas, aes(x = tfr_media, y = obs_por_millon_log)) +
  # Puntos con transparencia y tamaño fijo
  geom_point(aes(color = tfr_media), size = 3, alpha = 0.7) +
  # Línea de regresión lineal (OLS) para visualizar la asociación
  geom_smooth(method = "lm", se = TRUE, color = "#4682B4", linetype = "dashed") +
  
  # Etiquetado para identificar países clave
  geom_text_repel(
    aes(label = Etiqueta_pais), 
    box.padding = 0.5, 
    point.padding = 0.5,
    max.overlaps = Inf,
    segment.color = 'grey50',
    size = 3
  ) +
  
  scale_color_gradient(low = "#FFD700", high = "#B8860B", name = "TFR") +
  labs(
    title = "Frecuencia de Vigilancia (Normalizada) vs. Fertilidad Media",
    subtitle = "Log(Obs. de Aire/Millón de Hab.) en función de la TFR media",
    x = "Tasa de Fecundidad Total (TFR) Media",
    y = "Frecuencia de Vigilancia (Log de Obs. por Millón de Habitantes)"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold"))

Resultados del gráfico de dispersión:

El gráfico de dispersión muestra cómo se relacionan dos variables al colocar cada país como un punto en un plano.

Este gráfico evalúa la relación entre la Frecuencia de Reporte de Calidad del Aire (cuántas mediciones se envían, ajustado por población en el Eje Y) y la Tasa de Fecundidad Total (TFR) (en el Eje X).

Eje X: Cuanto más a la derecha está un país, mayor es su fertilidad media (TFR).
Eje Y: Cuanto más arriba está un país, mayor es su frecuencia de reporte de datos de aire por habitante (más datos envían).

Los puntos están distribuidos de forma caótica en todo el gráfico. No forman una línea o una curva clara que indique una relación.

La Línea de Tendencia (Regresión) que atraviesa la gráfica representa la tendencia media de la relación. Esta línea es casi perfectamente horizontal, lo que significa que el valor de la TFR (movernos a la derecha) no hace que la frecuencia de vigilancia suba ni baje (el punto Y se queda en el mismo nivel).Por lo tanto no existe una relación entre ambos.

Ejemplo: Hay países extremos, especialmente en la vigilancia, como España, Serbia, Eslovaquia y Alemania, que tienen niveles de vigilancia muy altos (puntos muy altos en el eje Y), pero sus TFRs son promedio o bajas.

Refutación de la Hipótesis Inicial: El gráfico sugiere que la hipótesis de que “una mayor vigilancia se asocia con niveles de fertilidad superiores” es probablemente falsa.

4.1.6 Test Estadístico: Regresión Lineal Robusta

Este apartado tiene como objetivo someter a prueba nuestra hipótesis central de que la frecuencia de la vigilancia ambiental se asocia con la fertilidad. Para ello, se emplea el Modelo de Mínimos Cuadrados Ordinarios (OLS).

Se utiliza el método de Errores Estándar Robustos (Huber-White, tipo HC1), que ajusta la varianza de los errores para obtener coeficientes y valores p más fiables, incluso bajo la presencia de heterocedasticidad.

La regresión se establece con la TFR media como variable dependiente ($Y$) y la variable Log(Obs. de Aire por Millón de Habitantes) ($\log(1+X_{Normalizada})$) como variable independiente o predictora.

# Filtro final del dataframe (Asegura que X e Y sean finitos)
datos_regresion <- df_vigilancia_tfr %>% 
    filter(is.finite(tfr_media), is.finite(obs_por_millon_log))

# Modelo OLS
modelo_final <- lm(tfr_media ~ obs_por_millon_log, data = datos_regresion)

resultados_robustez <- coeftest(
    modelo_final, 
    vcov = vcovHC(modelo_final, type = "HC1")
)

# Tabulación de resultados 
tabla_regresion_formato_simple <- tidy(resultados_robustez) %>%
  rename(
    T_Estadístico_Robusto = statistic,
    Valor_P_Robusto = p.value
  ) %>%
  mutate(
    # Redondeamos las estimaciones y errores a 4 decimales
    Estimación = round(estimate, 4),
    Error_Estandar_Robusto = round(std.error, 4),
    # Redondeamos el T-Estadístico a 3 decimales
    T_Estadístico_Robusto = round(T_Estadístico_Robusto, 3), 
    # Mantenemos el P-value con más precisión para su interpretación
    Valor_P_Robusto = round(Valor_P_Robusto, 6), 
    Significancia = case_when(
      Valor_P_Robusto < 0.001 ~ "Altamente significativo",
      Valor_P_Robusto < 0.01 ~ "Muy significativo",
      Valor_P_Robusto < 0.05 ~ "Significativo",
      TRUE ~ "No significativo"
    )
  ) %>%
  select(term, Estimación, Error_Estandar_Robusto, T_Estadístico_Robusto, Valor_P_Robusto, Significancia)

# Generación de la tabla
knitr::kable(
  tabla_regresion_formato_simple, 
  caption = "Resultados de Regresión: TFR Media ~ Log(Obs. de Aire por Millón de Hab.) con Errores Estándar Robustos (HC1)",
  # Alineación de las columnas (L=Left, R=Right)
  align = c('l', 'r', 'r', 'r', 'r', 'c') 
)

Resultados de Regresión: TFR Media ~ Log(Obs. de Aire por Millón de Hab.) con Errores Estándar Robustos (HC1)
term	Estimación	Error_Estandar_Robusto	T_Estadístico_Robusto	Valor_P_Robusto	Significancia
(Intercept)	1.5593	0.0312	49.949	0.00000	Altamente significativo
obs_por_millon_log	-0.0009	0.0029	-0.320	0.74951	No significativo

4.1.7 Conclusión del Objetico 1: Vigilancia vs Fertilidad

Los resultados de la regresión lineal robusta confirman lo visualizado en el análisis descriptivo. El coeficiente de la variable Frecuencia de Vigilancia normalizada (obs_por_millon_log) es $-0.0009$, lo que indica una relación prácticamente nula entre ambas variables. Más importante aún, el Valor P Robusto ($0.7495$) es significativamente superior al umbral de significancia de $0.05$. Por lo tanto, no existe evidencia estadística significativa para concluir que la intensidad de la vigilancia ambiental de la calidad del aire (ajustada por población) tenga algún efecto o asociación lineal con la Tasa de Fecundidad Total (TFR) entre los países de la muestra. La hipótesis inicial queda refutada en este primer objetivo.

4.2 Caracterizar la fertilidad media de los países europeos y comparación con los niveles de dióxido de nitrógeno

En esta segunda parte ampliamos el análisis pasando de la “vigilancia del aire” a la calidad del aire en sí, centrándonos en un contaminante clave: el dióxido de nitrógeno (NO₂). Este compuesto es un marcador típico de tráfico rodado y combustión urbana, y se ha vinculado con diversos efectos adversos para la salud, el objetivo aquí es evaluar la asociación agregada entre los niveles medios de NO₂ en países europeos y su Tasa de Fertilidad Total (TFR)..

Pregunta: ¿Existe una asociación entre los niveles medios de NO₂ en cada país europeo y su fertilidad media?

Hipótesis: Se espera una relación negativa, donde una mayor contaminación media por NO₂ se asocie con TFR más bajas, reflejando el entorno urbano e industrializado que suele caracterizar a los países con baja natalidad.

Como en el caso anterior, esta hipótesis es descriptiva, no causal; buscamos el signo de la asociación en los datos agregados.

4.2.1 Datos y construcción de variables

Para esta sección combinamos dos fuentes principales:

OMS — AAP 2022 (OMS_Datos.xlsx): Valores de NO₂ por ciudad.
Eurostat — Fertilidad (fert_resumen_nombres): TFR media por país.

# Leer directamente desde la hoja principal
oms <- read_excel("INPUT/DATA/OMS_Datos.xlsx", sheet = "AAP_2022_city_v9")

# Solo países de la región europea de la OMS
oms_europa <- oms %>%
  filter(`WHO Region` == "European Region")

# Variable 'pais' con nombres armonizados
oms_europa_limpio <- oms_europa %>%
  mutate(pais = recode_pais(`WHO Country Name`))

Tras armonizar y unir ambas tablas mediante left_join(), obtenemos un panel ciudad–país con valores de NO₂ y la fertilidad media asociada a cada país.
A partir de éste, construimos un resumen por país con las siguientes variables:

no2_media: media del NO₂ por país.
no2_mediana: mediana del NO₂.
n_ciudades: número de ciudades o localidades monitorizadas.
tfr_media: fertilidad media (Eurostat).

Este resumen nos permite trabajar con una unidad clara (el país) y explorar la relación entre contaminación y fertilidad.

# JOIN: tabla OMS ciudad-país + resumen de fertilidad por país
oms_europa_con_fert <- oms_europa_limpio %>%
  left_join(fert_resumen_nombres, by = "pais")

# Comprobar qué países de OMS no tienen fertilidad disponible
faltan_fert_oms <- oms_europa_limpio %>%
  distinct(pais) %>%
  anti_join(fert_resumen_nombres %>% distinct(pais), by = "pais")

# Resumen por país
pais_no2_fert <- oms_europa_con_fert %>%
  group_by(pais) %>%
  summarise(
    no2_media   = mean(`NO2 (μg/m3)`, na.rm = TRUE),
    no2_mediana = median(`NO2 (μg/m3)`, na.rm = TRUE),
    n_ciudades  = n_distinct(`City or Locality`),
    tfr_media   = first(tfr_media),   # viene del resumen de fertilidad (ya homog.)
    .groups = "drop"
  ) %>%
  filter(!is.na(no2_media), !is.na(tfr_media))

Del mismo modo, se recopilarán datos referentes a la clasificación por tipo de área de las estaciones de vigilancia de la calidad del aire, de modo que no solo tengamos información ordenada sobre los niveles de NO₂ en comparación con los de TFR, sino también de la procedencia de dicho gas.

4.2.2 Análisis Descriptivo I: Gráfico de Burbujas

Para visualizar la relación NO₂–fertilidad utilizamos un gráfico de burbujas:

Eje X: concentración media de NO₂ por país (no2_media).
Eje Y: fertilidad media (TFR, tfr_media).
Tamaño del punto: número de ciudades monitorizadas (n_ciudades).

La recta de regresión lineal (OLS, línea discontinua) se incluye para estimar visualmente la dirección y la fuerza de la asociación.

Este enfoque gráfico es fundamental para determinar si la hipótesis de una relación negativa entre NO2 y la TFR se manifiesta en los datos agregados europeos.

library(ggplot2)
library(ggrepel)

# Define el dataframe y las asignaciones a los ejes (Aesthetics)
ggplot(
  pais_no2_fert,
  aes(
    x    = no2_media,
    y    = tfr_media,
    size = n_ciudades
  )
) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, linetype = "dashed") +
  geom_text_repel(
    aes(label = pais),
    size        = 3,
    max.overlaps = 30
  ) +
  scale_size_continuous(name = "Nº de ciudades\nmonitorizadas") +
  labs(
    title    = "NO"[2]~" medio vs fertilidad media por país europeo",
    subtitle = "Tamaño del punto = nº de ciudades con datos de NO"[2]~"",
    x        = "NO"[2]~" medio (µg/m³, OMS)",
    y        = "TFR media (Eurostat)"
  ) +
  theme_minimal()

Resultados del gráfico de dispersión

En el grafico se observa:

Tendencia General (Línea Discontinua): La línea de regresión (la línea azul discontinua) es casi horizontal y muy plana. Esto significa que a medida que te mueves de izquierda a derecha (es decir, de países menos contaminados a países más contaminados por NO2), la TFR media (el Eje Y) apenas cambia.
Banda de Confianza (Sombra Gris): La banda de sombra gris alrededor de la línea es muy ancha. Esto indica que la tendencia observada es poco fiable.
Dispersión: Los puntos de los países están muy dispersos y no siguen un patrón claro a lo largo de la línea.

Se esperaba una relación negativa; es decir, que los países con mayores niveles de contaminación urbana por NO₂ tendrían una Tasa de Fecundidad Total (TFR) más baja.

El nivel medio de contaminación urbana por NO₂ en un país no tiene una asociación observable o significativa con su TFR media.

4.2.3 Análisis Descriptivo y Prueba Formal de la Asociación $\text{NO}_2$

Para cerrar la sección de NO₂ y dar una respuesta concluyente a la pregunta de investigación, se complementa el análisis visual con estadísticas descriptivas y una prueba de significancia de la correlación de Pearson.

Se presentan las estadísticas básicas de las variables clave del análisis (la concentración media de NO₂ por país y la TFR media) para entender su distribución y rango.

# Calcular estadísticas descriptivas para NO2 y TFR (usando los datos de la OMS)
no2_descriptivas <- pais_no2_fert %>%
  summarise(
    n_paises = n(),
    
    # Estadísticas de NO2 (Contaminación)
    no2_media_total = mean(no2_media, na.rm = TRUE),
    no2_sd = sd(no2_media, na.rm = TRUE),
    no2_min = min(no2_media, na.rm = TRUE),
    no2_max = max(no2_media, na.rm = TRUE),
    
    # Estadísticas de TFR (Fertilidad)
    tfr_media_total = mean(tfr_media, na.rm = TRUE),
    tfr_sd = sd(tfr_media, na.rm = TRUE),
    tfr_min = min(tfr_media, na.rm = TRUE),
    tfr_max = max(tfr_media, na.rm = TRUE)
  )

Resultados de la Prueba de Asociación

El análisis descriptivo de las variables en los 36 países revela grandes diferencias en la variación:

NO₂ Media (Contaminación): La media es de $18.73\ \mu \text{g}/\text{m}^3$. Existe una amplia variación en los niveles de contaminación (rango de $6.77$ a $29.98\ \mu \text{g}/\text{m}^3$), lo que significa que los países de la muestra tienen niveles de NO₂ muy distintos.
TFR Media (Fertilidad): La media es de $1.55$ nacimientos por mujer. La variación es mucho más baja (rango de $1.24$ a $1.96$), confirmando que los patrones de baja natalidad son la norma homogénea en Europa.

En resumen: La contaminación por NO₂ varía mucho más que la TFR. El análisis de correlación determinará si esta gran variación en NO₂ tiene algún poder para explicar la pequeña variación en la fertilidad.

4.2.4 Prueba de Significancia de la Correlación

Para determinar si la correlación débil observada en el gráfico es estadísticamente significativa o simplemente un producto del azar en esta muestra, se utiliza una prueba de significancia de la correlación de Pearson (cor.test).

# Realizar la prueba de significancia de la correlación de Pearson
prueba_correlacion <- cor.test(pais_no2_fert$no2_media, pais_no2_fert$tfr_media, method = "pearson")

4.2.5 Conclusión del Objetico 2: NO₂ vs Fertilidad

El análisis de correlación formal entre los niveles de NO₂ medio y la Tasa de Fecundidad Total (TFR) en 36 países europeos revela una ausencia de asociación. El coeficiente de Pearson ($\rho$) de $0.0724$ indica una relación lineal prácticamente nula, que es contraria a la hipótesis inicial de una relación negativa. Además, esta correlación resulta ser estadísticamente no significativa, con un Valor P de $0.675$ (muy superior al umbral de $0.05$).

En consecuencia, se rechaza la hipótesis de que la contaminación atmosférica urbana (medida por el NO₂) es un factor determinante de la fertilidad a nivel agregado. El NO₂ no actúa como un predictor sistemático de la TFR; en su lugar, las variaciones en la natalidad están impulsadas por factores macroeconómicos, culturales y de políticas públicas que no están correlacionados con los niveles medios de este contaminante.

4.3 Descubrir si la relación entre el dióxido de carbono y la fertilidad es directamente proporcional

En este tercer objetivo analizamos el papel de un contaminante global: las Emisiones de $\text{CO}_2$ per cápita.

Exploramos la asociación entre las emisiones de CO₂ y la Tasa de Fecundidad Total (TFR) para evaluar si las estructuras económicas de alto impacto ambiental se correlacionan con los patrones de baja natalidad.

Pregunta: ¿Existe una relación entre las Emisiones de CO₂ per cápita en los países europeos y su fertilidad media?

Hipótesis: Se espera una relación negativa. Los países con mayores emisiones de CO₂ (altamente industrializados/desarrollados) presentarán TFR bajas.

4.3.1 Datos y construcción de variables

Partimos de 2 ficheros:

JSON (CO2.json) que contiene las emisiones anuales de $\text{CO}_2$ (en kilotoneladas) por país y por año.
Eurostat (fert_resumen_nombres) para la TFR.

# Leemos el JSON como Dataframe, no se por qué pero daba errores sino
dioxidocarbono_df <- jsonlite::fromJSON(
  "INPUT/DATA/CO2.json",
  simplifyDataFrame = TRUE  # lo aplana a data.frame directamente
) %>%
  as_tibble() %>%
  clean_names()

# Tabla co2_limpio (país-año)
co2_limpio <- dioxidocarbono_df %>%
  transmute(
    pais = recode_pais(country_name),
    anio = as.integer(year),
    co2_kt = as.numeric(co2_emissions_kt)
  ) %>%
  group_by(pais, anio) %>%
  summarise(
    co2_kt = sum(co2_kt, na.rm = TRUE),
    .groups = "drop"
  )
# Panel CO2 + FERTILIDAD (país–año)
co2_fert_panel <- co2_limpio %>%
  inner_join(fert_paises, by = c("pais", "anio"))

# Resumen por país
co2_resumen_pais <- co2_limpio %>%
  group_by(pais) %>%
  summarise(
    n_anios_co2  = n_distinct(anio),
    anio_min_co2 = min(anio, na.rm = TRUE),
    anio_max_co2 = max(anio, na.rm = TRUE),
    co2_kt_media = mean(co2_kt, na.rm = TRUE),
    co2_kt_total = sum(co2_kt, na.rm = TRUE),
    .groups = "drop"
  )

co2_fert_resumen_pais <- co2_resumen_pais %>%
  full_join(fert_resumen_nombres, by = "pais") %>%
  arrange(pais)

4.3.2 Análisis Descriptivo I: Gráfico de ranking tipo “lollipop”

La primera estrategia visual para abordar la relación CO₂-Fertilidad es un gráfico de ranking que ordena a los países por sus Emisiones Medias de CO₂ a lo largo del periodo de estudio.

Este formato de visualización, permite:

Establecer un ranking claro para identificar rápidamente a los mayores y menores emisores CO₂ total.
Superponer la TFR: La Tasa de Fecundidad Total (tfr_media) se representa mediante el color y el tamaño del punto al final de cada barra.

Al ordenar los países por CO₂ y colorearlos por TFR, se puede evaluar visualmente si existe una agrupación sistemática. Si la hipótesis de una relación negativa fuera cierta, esperaríamos que los países con las barras de CO₂ más largas tuvieran puntos más pequeños y de color menos intenso (TFR baja), y viceversa.

co2_lolli <- co2_fert_resumen_pais %>%
  filter(!is.na(co2_kt_media),
         !is.na(tfr_media))

ggplot(
  co2_lolli,
  aes(
    x = co2_kt_media,
    y = fct_reorder(pais, co2_kt_media)
  )
) +
  # palito
  geom_segment(
    aes(
      x    = 0,
      xend = co2_kt_media,
      y    = fct_reorder(pais, co2_kt_media),
      yend = fct_reorder(pais, co2_kt_media)
    ),
    linewidth = 0.6,
    alpha     = 0.6
  ) +
  # punto, coloreado por TFR
  geom_point(
    aes(
      color = tfr_media,
      size  = tfr_media
    ),
    alpha = 0.9
  ) +
  scale_size_continuous(name = "TFR media") +
  scale_color_viridis_c(name = "TFR media", option = "C") +
  labs(
    title    = "Ranking de emisiones medias de CO2 y fertilidad en Europa",
    subtitle = "Línea = emisiones medias de CO2 (kt) | Color/tamaño = TFR media",
    x        = "CO2 medio (kt, total país)",
    y        = NULL
  ) +
  theme_minimal()

Resultados del Gráfico de Lollipop

El análisis del ranking de CO₂ total frente a la TFR media no respalda la hipótesis de una relación negativa (más CO₂ significa menos fertilidad). El gráfico muestra que los mayores emisores (ej., Polonia) tienen niveles de TFR ($\approx 1.5$) que son similares a la media de la muestra, y los menores emisores (ej., Malta) pueden registrar las TFR más bajas, rompiendo cualquier patrón sistemático. Por lo tanto, el volumen total de emisiones de CO₂ de un país no es un factor que diferencie significativamente los patrones de natalidad en Europa.

4.3.3 Análisis Descriptivo II: Gráfico animado CO₂–TFR por año

Para un análisis más dinámico, se presenta la relación CO₂-Fertilidad en un diagrama de dispersión animado, donde cada fotograma representa un año concreto del periodo de estudio.

La animación permite:

Observar la Tendencia Temporal: Cada fotograma muestra el diagrama de dispersión entre las emisiones totales de CO₂ (co2_kt) y la TFR (tfr) para un año específico.
Evaluar la Fuerza Anual: Se ajusta una recta de regresión lineal con su banda de confianza a los datos de cada año. La inspección visual de si la recta es plana o inclinada permite ver la tendencia dominante en ese momento.
Cuantificación Métrica: En una esquina de la animación, se añaden los coeficientes estimados (la pendiente “Beta” y el $\text{R}^2$), junto con el número de países, permitiendo cuantificar la fuerza de la relación en ese año específico.

De esta forma, la animación permite ver si la relación cambia con el tiempo y si, en algún periodo concreto, la hipótesis de una tendencia clara (positiva o negativa) llega a ser temporalmente significativa.

co2_fert_anim <- co2_fert_panel %>%
  filter(!is.na(co2_kt),
         !is.na(tfr))

reg_stats <- co2_fert_anim %>%
  group_by(anio) %>%
  summarise(
    n    = n(),
    beta = coef(lm(tfr ~ co2_kt))[2],
    r2   = summary(lm(tfr ~ co2_kt))$r.squared,
    .groups = "drop"
  )

library(gganimate)

p_anim <- ggplot(
  co2_fert_anim,
  aes(
    x = co2_kt,
    y = tfr
  )
) +
  # puntos por país
  geom_point(
    aes(color = pais),
    alpha = 0.7,
    size  = 2,
    show.legend = FALSE
  ) +
  # nombres de algunos países
  ggrepel::geom_text_repel(
    data = ~ dplyr::filter(
      .x,
      co2_kt == max(co2_kt) | tfr == max(tfr)
    ),
    aes(label = pais),
    size         = 3,
    max.overlaps = 30,
    show.legend  = FALSE
  ) +
  # recta de regresión por año
  geom_smooth(
    method  = "lm",
    se      = TRUE,
    linetype = "dashed",
    color   = "deeppink"
  ) +
  # texto con stats por año (beta, R², n)
  geom_text(
    data = reg_stats,
    aes(
      x     = -Inf,
      y     = Inf,
      label = sprintf("Beta = %.3f   R2 = %.2f   n = %d", beta, r2, n)
    ),
    hjust        = -0.05,
    vjust        = 1.2,
    size         = 3.2,
    inherit.aes  = FALSE
  ) +
  scale_x_continuous(labels = scales::comma) +
  labs(
    title    = "Relación entre emisiones de CO₂ y fertilidad en Europa",
    subtitle = "Año: {frame_time}",
    x        = "CO2 total del país (kt)",
    y        = "Tasa global de fecundidad (TFR)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.minor = element_blank(),
    plot.title       = element_text(face = "bold")
  ) +
  transition_time(anio) +
  ease_aes("linear")

# Ver la animación en el viewer
anim <- animate(
  p_anim,
  nframes  = 150,
  fps      = 10,
  renderer = gifski_renderer()
)

anim

Resultados del gráfico animado

El gráfico de dispersión, que compara las CO₂ medias frente a la TFR medias, demuestra que la asociación es nula y la hipótesis se rechaza visualmente.

Nula Tendencia: La nube de puntos está extremadamente dispersa. La línea de regresión es casi plana (horizontal), lo que indica que el CO₂ no tiene capacidad para explicar la variación en la TFR.
Hipótesis Refutada: La hipótesis de una relación negativa (“más CO₂ = menos TFR”) queda descartada. Los países se distribuyen de forma aleatoria, confirmando que las variables socioeconómicas asociadas al desarrollo (que el CO₂ representa) no determinan la fertilidad de forma sistemática.

4.3.4 Análisis Descriptivo III: Gráfico de Dispersión

Tras analizar el ranking y la evolución temporal, la siguiente estrategia visual es un Gráfico de Dispersión Agregado.

Este formato de visualización, que compara las CO₂ medias frente a la TFR medias, permite:

Visualizar la Nube de Puntos: Muestra la distribución de todos los países en el plano CO₂-TFR, revelando de inmediato si la asociación es dispersa (débil) o concentrada (fuerte).
Evaluar la Tendencia: La recta de regresión ajustada y su banda de confianza indican la dirección (positiva o negativa) y la robustez de la tendencia lineal general.
Contrastar la Hipótesis: La hipótesis de una relación negativa requiere que la recta de regresión tenga una pendiente descendente. Si la recta es plana, la hipótesis no se verifica.

Al inspeccionar la pendiente de la recta y la dispersión de los puntos, se obtiene la última confirmación visual antes de pasar a cuantificar la asociación con el Valor P.

# Calculamos la correlación agregada para mostrarla en el gráfico
co2_correlacion_agregada <- co2_fert_resumen_pais %>%
  filter(!is.na(co2_kt_media), !is.na(tfr_media)) %>%
  summarise(
    cor_pearson = cor(co2_kt_media, tfr_media, use = "complete.obs")
  )

# Creamos el gráfico de dispersión estático
ggplot(
  co2_fert_resumen_pais %>% filter(!is.na(co2_kt_media), !is.na(tfr_media)),
  aes(
    x = co2_kt_media,
    y = tfr_media,
    label = pais # Se añaden etiquetas de país para identificar outliers
  )
) +
  geom_point(alpha = 0.7, color = "#9B59B6", size = 2) + # Puntos de datos
  geom_smooth(method = "lm", se = TRUE, linetype = "dashed", color = "#3498DB") + # Línea de regresión
  ggrepel::geom_text_repel(
    size = 3,
    max.overlaps = 30
  ) +
  labs(
    title = "CO"[2]~" medio total vs. Fertilidad media en Europa (Análisis Agregado)",
    subtitle = paste0("Correlación de Pearson (ρ): ", round(co2_correlacion_agregada$cor_pearson, 3)),
    x = "CO"[2]~"  medio por país (Kilotoneladas, kt)",
    y = "TFR media (Eurostat)"
  ) +
  theme_minimal()

Resultados del Gráfico de Dispersión

El gráfico de dispersión confirma que la relación CO₂-Fertilidad es prácticamente nula:

Correlación ($\rho$) Nula: El coeficiente de correlación de $\mathbf{-0.043}$ está pegado a cero, lo que significa que el CO₂ no tiene poder de explicación.
Línea Plana: La línea de tendencia es casi perfectamente horizontal, indicando que la TFR no cambia sistemáticamente, sin importar si el país tiene mucho o poco CO₂.

La hipótesis inicial de una relación negativa (a más CO₂, menos fertilidad), por lo que se rechaza visualmente.

4.3.5 Prueba de significancia

El paso final es aplicar la Prueba de Correlación de Pearson (cor.test) para determinar:

Cuantificación Final: El valor exacto del coeficiente de correlación ($\rho$) para todo el periodo.
Significancia Estadística: El Valor P. Si este valor es superior a $0.05$, se concluye que la correlación observada es un producto del azar y no es válida para la población.

Este análisis proporcionará el veredicto estadístico final sobre la hipótesis del Objetivo 3.

# Aseguramos que solo incluimos países con ambos datos (CO2 medio y TFR media)
datos_cor_co2 <- co2_fert_resumen_pais %>%
  filter(!is.na(co2_kt_media), !is.na(tfr_media))

# Realizar la prueba de significancia de la correlación de Pearson
prueba_cor_co2 <- cor.test(datos_cor_co2$co2_kt_media, datos_cor_co2$tfr_media, method = "pearson")

4.3.6 Conclusiones del Objetivo 3: CO₂ vs fertilidad

El análisis formal de la asociación entre las Emisiones medias de CO₂ y la Tasa de Fecundidad Total (TFR) en los países europeos confirma las impresiones visuales obtenidas en los gráficos. El coeficiente de Pearson ($\rho$) es de $\mathbf{-0.043}$, un valor prácticamente nulo que indica que las emisiones de CO₂ no tienen una relación lineal con la TFR. Este resultado contradice la hipótesis de una relación negativa fuerte. La falta de relación se confirma mediante la prueba de significancia: el Valor P es de $\mathbf{0.8839}$, siendo superior al umbral de $\mathbf{0.05}$. Por lo tanto, se concluye formalmente que no existe una asociación estadísticamente significativa entre el CO₂ total y la TFR media en Europa.

La hipótesis de una relación entre el desarrollo industrial y la baja natalidad, medida a través de las emisiones de CO₂, se rechaza. Las variaciones en la TFR están impulsadas por otros factores socioeconómicos y culturales que el CO₂ no logra capturar de manera sistemática; además, sabiendo que los niveles de este gas son directamente proporcionales a los de tráfico, es evidente que una reducción de éste redduciría considerablemente este factor contaminante.

Por tanto, el tercer objetivo confirma la idea general de los anteriores:
cuando miramos a nivel país, con indicadores agregados, la relación entre contaminación (o emisiones) y fertilidad en Europa parece ser, como mínimo, muy débil y fuertemente confundida por otras dimensiones estructurales.

4.4 Caracterizar la importancia del PM₂.₅ dentro de la vigilancia del aire y su relación con la fertilidad

En este último objetivo profundizamos en un contaminante concreto, el material particulado fino (PM₂.₅), que se ha asociado con efectos adversos sobre la salud reproductiva, el embarazo y la mortalidad perinatal. A diferencia de los objetivos anteriores, centrados en la intensidad global de la vigilancia o en contaminantes específicos como el NO₂ o el CO₂, aquí nos interesa el peso específico de las particulas de PM₂.₅ dentro la monitorización de la calidad del aire en cada país. Este indicador mide el tamaño de las partículas en el aire, tomando en cuenta aquellas de 2,5 micrómetros o menos; cuanto más pequeñas son estas partículas, más profundo llegan a los pulmones. Es claro que personas con patologías respiratorias sufrirán más la presencia de un PM₂.₅ alto, aunque a largo plazo toda la población se acaba viendo afectada por ellas. Por ello, el PM₂.₅ está regulado por la Agencia Europea de Medio ambiente (EEA) y por la OMS, que recomiendan niveles bajo de exposición por su impacto demostrado en la salud.

La lógica es la siguiente: Un país que dedica una mayor proporción de sus observaciones al PM₂.₅ estaría, en principio, más orientado a monitorizar riesgos finos para la salud, lo que podría estar asociado indirectamente con entornos más saludables y una planificación familiar más favorable.

Pregunta: ¿Existe una relación entre el peso relativo del PM₂.₅ en la vigilancia del aire y la fertilidad media (TFR) de los países europeos?

Hipótesis: Si el PM₂.₅ es un marcador relevante de riesgos ambientales percibidos o efectivos, cabría esperar que los países que dedican una mayor porción de sus observaciones a este contaminante presenten entornos más saludables o más vigilados, lo que podría asociarse con una fertilidad algo más elevada. De forma descriptiva, planteamos la hipótesis de que una mayor proporción de observaciones de PM₂.₅ se asocia con una TFR media más alta, aunque sin asumir causalidad.

4.4.1 Datos y construcción de variables

En este objetivo analizamos el lugar que ocupa el contaminante PM₂.₅ dentro de la monitorización de calidad del aire en Europa.

Partimos del fichero de registros de calidad del aire (aire), que contiene observaciones por país, año y tipo de contaminante. A partir de este conjunto:

Se armonizan los nombres de país (pais), al igual que en los objetivos anteriores.
Se convierten los años a formato numérico (anio).
Se crea la variable lógica es_pm25, que toma valor TRUE cuando el contaminante corresponde a PM₂.₅ (incluyendo variantes como "PM2.5", "PM25", "PM2p5").

Con esta base construimos el resumen por país pm25_share_pais, que contiene:

n_obs_total: número total de observaciones,
n_obs_pm25: observaciones PM₂.₅,
prop_pm25: proporción de observaciones dedicadas a PM₂.₅,
n_anios: número de años con observaciones.

Después unimos este resumen con fert_resumen_nombres para obtener pm25_share_tfr, que contiene tanto la proporción de PM₂.₅ como la fertilidad media (tfr_media) por país.

Finalmente, unimos este resumen con fert_resumen_nombres para obtener pm25_share_tfr. Se aplica un filtrado final para asegurar que solo se incluyan los países que poseen datos válidos tanto de la proporción de PM₂.₅ como de la Tasa de Fecundidad Total (TFR), garantizando la integridad de la muestra para el análisis.

# Metadatos de estaciones EEA (dataframe interno)

stations_eea <- EEAaq_get_dataframe("stations") %>%
  clean_names() %>%
  mutate(
    pais_raw = str_trim(country),
    pais = recode_pais(country)
  )

# Renombramos columna problemática
stations_meta <- stations_eea %>%
  rename(air_quality_station_eoi_code = air_quality_station_eo_i_code) %>%
  select(
    air_quality_station_eoi_code,
    pais,
    air_quality_station_area,
    air_quality_station_type
  )

# Cálculo de proporción de PM2.5 por país
pm25_share_pais <- aire %>%
  mutate(
    pais = recode_pais(country),
    anio = as.integer(year),
    es_pm25 = air_pollutant %in% c("PM2.5", "PM25", "PM2p5")
  ) %>%
  group_by(pais) %>%
  summarise(
    n_obs_total = n(),
    n_obs_pm25  = sum(es_pm25, na.rm = TRUE),
    prop_pm25   = n_obs_pm25 / n_obs_total,
    n_anios     = n_distinct(anio),
    .groups     = "drop"
  ) %>%
  arrange(desc(prop_pm25))


# Unión PM2.5 + TFR
pm25_share_tfr <- pm25_share_pais %>%
  left_join(fert_resumen_nombres, by = "pais") %>%
  filter(!is.na(tfr_media), !is.na(prop_pm25))

No hay que olvidar que el análisis inicial de la proporción total de PM₂.₅ (prop_pm25) puede ser engañoso, ya que una medición baja puede deberse a que el país prioriza otras variables como las analizadas anteriormente, o bien tiene estaciones de vigilancia concentradas en zonas de tráfico o industriales.

4.4.2 Análisis Descriptivo I: Gráfico de Barras Ordenadas por Proporción de $\text{PM}_{2.5}$

Para examinar de manera descriptiva la posible asociación entre la prioridad técnica que un país otorga al PM₂.₅ y su fertilidad, se ha optado por un Gráfico de Barras Ordenadas.

Evaluación de Prioridad: La altura de la barra (Eje Y) representa el porcentaje de todas las observaciones de vigilancia de calidad del aire que un país dedica al PM₂.₅. Los países se ordenan progresivamente en el Eje X de menor a mayor prioridad.
Relación con la TFR: El color de la barra codifica la Tasa de Fecundidad Total (TFR) media.

Al ordenar la variable independiente ($\text{Prop}_{\text{PM}_{2.5}}$) en el eje, el gráfico permite identificar si existe un patrón visual donde el color se intensifica (o se aclara) a medida que las barras se hacen más altas. Si la hipótesis fuera correcta, observaríamos un claro gradiente de color que va de claro a oscuro (mayor TFR) a lo largo del eje X. Este análisis descriptivo sirve como paso previo indispensable antes de confirmar la relación mediante la regresión robusta.

# Ordenamos países por prop_pm25 y preparamos factor
pm25_circ <- pm25_share_tfr %>%
  arrange(prop_pm25) %>% # Ordenamos ascendentemente para el eje X
  mutate(
    pais      = factor(pais, levels = pais),
    label_pm  = scales::percent(prop_pm25, accuracy = 1)
  )

ggplot(
  pm25_circ,
  aes(
    x    = pais,
    y    = prop_pm25,
    fill = tfr_media
  )
) +
  # Barras (geom_col), coloreadas por la TFR media
  geom_col(aes(fill = tfr_media), stat = "identity") +
  
  # Escala de color para representar la fertilidad (misma escala que Obj 1)
  scale_fill_gradient(
    low = "peachpuff", 
    high = "red", 
    space = "Lab", 
    name = "TFR media"
  ) +
  
  # Etiquetas y formato del eje Y
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    expand = expansion(mult = c(0, .1))
  ) +
  
  labs(
    title = expression("Asociación entre Prioridad de Vigilancia PM"[2.5]~"y Fertilidad"),
    subtitle = expression("Países ordenados por Proporción de Observaciones PM"[2.5]~"(Eje Y) | Color = TFR media"), 
    x = NULL, 
    y = expression("Proporción de Observaciones PM"[2.5]~"sobre el Total"),
    fill = "TFR media"
  ) +
  
  theme_classic() +
  theme(
    axis.text.x = element_text(angle = 60, vjust = 1, hjust = 1, size = 8),
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

Resultado del gráfico de barras

El análisis del Gráfico de Barras Ordenadas revela dos mensajes clave sobre la prioridad de monitorización del PM₂.₅ en el contexto de la fertilidad europea:

Baja Prioridad Media: En la mayoría de los países europeos, el esfuerzo de vigilancia dedicado al PM₂.₅ es minoritario, observando una mayor proporción por debajo del 10% del total. Esta heterogeneidad va desde valores marginales (alrededor del 2% en Luxemburgo) hasta la prioridad más alta (cercana al 18% en Noruega).
Ausencia de Patrón Monotónico: Crucialmente, al examinar el color (TFR media) a lo largo del eje X (prioridad de PM₂.₅), no se observa una tendencia visual coherente que apoye la hipótesis inicial. Los tonos más oscuros (TFR alta) no se concentran de forma sistemática en el extremo derecho (mayor proporción de PM₂.₅). Países como Irlanda o Noruega, que tienen una alta prioridad al PM₂.₅, muestran TFR medias (tonos intermedios), mientras que países con TFRs relativamente altas aparecen dispersos a lo largo de todo el rango de prioridad de vigilancia.

En conjunto, la evidencia descriptiva no sugiere la existencia de una asociación positiva entre el peso del PM₂.₅ en la vigilancia de la calidad del aire y la TFR media. Si existiera un efecto, se vería un cambio gradual y ordenado en el color. Como no lo vemos, el gráfico nos dice de forma preliminar que nuestra hipótesis no se sostiene.

4.4.3 Análisis Descriptivo II: Gráfico de Dispersión

El gráfico de dispersión permite evaluar si existe algún patrón funcional entre la $\text{Prop}_{\text{PM}_{2.5}}$ y la TFR, además de identificar posibles países atípicos (outliers).

Variable Predictora ($X$): Proporción de Observaciones PM₂.₅ (prop_pm25).

Variable Dependiente ($Y$): Tasa de Fecundidad Total Media (tfr_media).

ggplot(pm25_share_tfr, aes(x = prop_pm25, y = tfr_media)) +

  geom_point(aes(color = tfr_media), size = 3) +

  # Añadimos una línea de tendencia (regresión lineal) para guiar la vista

  geom_smooth(method = "lm", se = FALSE, color = "darkgrey", linetype = "dashed") +

  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +

  labs(

    title = "Relación entre la Proporción de Vigilancia PM"[2.5]~" y la Fertilidad Media",

    subtitle = "Cada punto es un país, la línea discontinua es la tendencia media (OLS)",

    x = "Proporción de Observaciones dedicadas a PM"[2.5]~" (Prop_PM"[2.5]~")",

    y = "Tasa de Fecundidad Total Media (TFR)",

    color = "TFR media"

  ) +

  theme_minimal() +

  theme(plot.title = element_text(face = "bold", hjust = 0.5))

Resultado del gráfico de dispersión

Nuestra hipótesis era que, si el PM₂.₅ es un riesgo importante y bien vigilado, los países con más prioridad de PM₂.₅ (puntos a la derecha) tendrían una TFR más alta (puntos más arriba y de color más claro/azul claro).Es decir, deberíamos ver una línea de tendencia claramente inclinada hacia arriba.

Sin embargo, lo que observamos es: - Los puntos están esparcidos por todo el gráfico, como si no tuvieran una dirección clara. Hay países con baja prioridad (PM₂.₅) que tienen TFR alta, y países con alta prioridad (PM₂.₅) que tienen TFR media.

La línea discontinua (que muestra la tendencia promedio) está casi completamente plana. Es cierto que tiene una ligera inclinación ascendente, pero es insignificante.
Los países europeos están casi todos agrupados en el rango de TFR entre 1.4 y 1.7 (tonos azules oscuros e intermedios). La prioridad del PM₂.₅ no parece sacarlos de ese rango.

En conclusión, no existe una relación visible entre el esfuerzo que un país dedica a medir el PM₂.₅ y su Tasa de Fecundidad Total (TFR). La línea de tendencia es tan plana que nos dice que, sin importar cuánto se preocupe un país por medir este contaminante (eje X), su fertilidad (eje Y) se mantiene casi igual.

4.4.4 Test Estadístico: Regresión Lineal Robusta

Tras el análisis descriptivo visual (gráfico de barras y de dispersión), que sugirió la ausencia de una asociación coherente entre la prioridad de vigilancia del PM₂.₅ ($\text{Prop}_{\text{PM}_{2.5}}$) y la fertilidad media (TFR), es necesario aplicar un test estadístico que confirme o rechace formalmente la hipótesis nula ($\text{H}_0$: no hay relación).

Para ello, utilizamos el modelo de Regresión Lineal Robusta (OLS con Errores Estándar Huber-White). Este método es el más adecuado para datos de corte transversal por país, ya que permite estimar el coeficiente de asociación ($\beta_1$) y su significancia ($p$-valor), al tiempo que neutraliza el riesgo de heterocedasticidad (varianza no constante de los errores) que es frecuente en este tipo de datos agregados.

El modelo nos permitirá determinar si el impacto de la $\text{Prop}_{\text{PM}_{2.5}}$ sobre la TFR es estadísticamente diferente de cero. Si el $p$-valor es alto ($p > 0.05$), concluiremos que la asociación es tan débil que no es fiable.

# Filtro final del dataframe
datos_regresion_pm25 <- pm25_share_tfr %>%
    filter(is.finite(tfr_media), is.finite(prop_pm25))

# Modelo OLS: TFR media ~ Proporción de PM₂.₅
modelo_pm25 <- lm(tfr_media ~ prop_pm25, data = datos_regresion_pm25)

# Aplicar Errores Estándar Robustos (Huber-White, tipo HC1)
resultados_pm25_robustez <- coeftest(
    modelo_pm25,
    vcov = vcovHC(modelo_pm25, type = "HC1")
)

# Tabulación de resultados con formato de texto 
tabla_regresion_pm25 <- tidy(resultados_pm25_robustez) %>%
  rename(
    T_Estadístico_Robusto = statistic,
    Valor_P_Robusto = p.value
  ) %>%
  mutate(
    Estimación = round(estimate, 4),
    Error_Estandar_Robusto = round(std.error, 4),
    T_Estadístico_Robusto = round(T_Estadístico_Robusto, 3), 
    Valor_P_Robusto = round(Valor_P_Robusto, 6), 
    Significancia = case_when(
      Valor_P_Robusto < 0.001 ~ "Altamente significativo",
      Valor_P_Robusto < 0.01 ~ "Muy significativo",
      Valor_P_Robusto < 0.05 ~ "Significativo",
      TRUE ~ "No significativo"
    )
  ) %>%
  select(term, Estimación, Error_Estandar_Robusto, T_Estadístico_Robusto, Valor_P_Robusto, Significancia)

# Generación de la tabla
knitr::kable(
  tabla_regresion_pm25,
  caption = "Resultados de Regresión: TFR Media ~ Proporción de Vigilancia PM₂.₅ con Errores Estándar Robustos (HC1)",
  align = c('l', 'r', 'r', 'r', 'r', 'l')
)

Resultados de Regresión: TFR Media ~ Proporción de Vigilancia PM₂.₅ con Errores Estándar Robustos (HC1)
term	Estimación	Error_Estandar_Robusto	T_Estadístico_Robusto	Valor_P_Robusto	Significancia
(Intercept)	1.4573	0.0612	23.817	0.0000	Altamente significativo
prop_pm25	1.4253	0.8398	1.697	0.0983	No significativo

4.4.5 Conclusión del Objetivo 4: $\text{PM}_{2.5}$ en la vigilancia vs fertilidad

El análisis descriptivo, tanto en el gráfico de barras como en el de dispersión , demostró visualmente que la Tasa de Fecundidad (TFR) no sigue un patrón claro a medida que aumenta la prioridad de vigilancia del contaminante $PM₂.₅ en un país. La línea de tendencia en el gráfico de dispersión se mostró prácticamente plana, sugiriendo que la relación era débil o inexistente. Al llevar a cabo el Test Estadístico de Regresión Lineal Robusta, esta observación se confirmó formalmente. A pesar de que la estimación mostró una ligera inclinación positiva (apoyando superficialmente la hipótesis inicial), el Valor P fue de $\text{0.0983}$, el cual es superior al umbral de confianza del $\text{5\%}$ ($p > 0.05$).

En términos sencillos, el resultado nos indica que, aunque la relación que medimos apunta en la dirección de nuestra hipótesis, es estadísticamente insignificante; no es lo suficientemente fuerte para descartar que sea pura casualidad. Por lo tanto, la hipótesis de que una mayor prioridad en la monitorización de PM₂.₅ se asocia a una mayor fertilidad queda rechazada. Esto concluye que la variabilidad en la TFR entre países europeos está impulsada por factores socioeconómicos y demográficos de mayor peso, y no por las decisiones técnicas relativas a la estrategia de vigilancia de la calidad del aire.

5 Conclusiones

El objetivo principal de este proyecto era ver si la contaminación ambiental (NO₂ y PM₂.₅) y las emisiones económicas (CO₂) tienen alguna relación con cuántos hijos nacen en Europa (la Tasa de Fecundidad Total o TFR). Los resultados, después de analizar los tres indicadores, son contundentes: no encontramos ninguna relación lineal estadísticamente significativa entre estas variables a nivel de país.

Para el contaminante urbano NO₂, la correlación ($\rho$) fue casi cero ($\approx 0.07$) y, lo más importante, el Valor P de $\mathbf{0.675}$ nos dice que esta pequeña tendencia es fruto del azar.
Para las emisiones de CO₂ (que representan el desarrollo industrial de cada país), la correlación fue igualmente nula ($\rho \approx -0.043$) y el Valor P de $\mathbf{0.8839}$ confirma que no hay conexión real.

En resumen, la hipótesis de que un país más contaminado o más industrializado automáticamente tiene menos bebés se rechaza en este estudio.

La clave de este resultado radica en la naturaleza de la TFR. El estudio demostró que la TFR media es muy estable y se mueve en un rango muy estrecho (alrededor de $1.55$ nacimientos por mujer) en casi toda Europa. Esta estabilidad contrasta enormemente con la gran diferencia que existe en los niveles de NO₂ y CO₂ entre un país y otro. Si la TFR variase mucho, quizás el CO₂ o el NO₂ podrían influir, pero como casi todos los países están en la misma situación demográfica, la contaminación no tiene espacio para ser un factor diferenciador. Esto confirma que los patrones de natalidad están dominados por variables mucho más poderosas que la contaminación atmosférica media.

Aunque las conclusiones son claras a nivel agregado, es vital mencionar la limitación principal de nuestro análisis:

Escala Agregada: Todos los indicadores de contaminación y TFR se promediaron a nivel de país. Esta metodología es ideal para ver tendencias macroeconómicas, pero es muy mala para capturar efectos biológicos directos.

Para futuros análisis, sería necesario utilizar datos con mucha más resolución, comparando la exposición directa a contaminantes (ej., el PM₂.₅ en el código postal) con las tasas de natalidad o los datos de salud a nivel micro (municipio o región).

6 Referencias

Eurostat. (s. f.). Fertility statistics (TFR by country and year). Oficina de Estadística de la Unión Europea.

Air Quality e-Reporting (AQ e-Reporting). (2022, August 5). https://www.eea.europa.eu/en/datahub/datahubitem-view/3b390c9c-f321-490a-b25a-ae93b2ed80c1

Emisiones globales de CO2. (2025, August 18). https://data.unesco.org/explore/dataset/co2001/information/?flg=es-es

World Health Organization. (2022). Ambient Air Pollution Database (AAP 2022). WHO, Department of Environment, Climate Change and Health.

Carré, J., Gatimel, N., Moreau, J., Parinaud, J., & Léandri, R. (2017). Does air pollution play a role in infertility? A systematic review. Environmental Health, 16(1), 82.

Landrigan, P. J., Fuller, R., Acosta, N. J. R., et al. (2018). The Lancet Commission on pollution and health. The Lancet, 391(10119), 462–512.

WHO. (2021). WHO global air quality guidelines: Particulate matter (PM₂.₅ and PM₁₀), ozone, nitrogen dioxide, sulfur dioxide and carbon monoxide. World Health Organization.

Impacto de la vigilancia y calidad del aire sobre la fertilidad en la población europea

Daniel Marina de la Cal, Patricia Ibarrondo Revilla y Nerea Yi Estepar Calvo

2025-11-24